clear all
nx=2;
ny=2;%even

a1=[1;0];
a2=[1/2;sqrt(3)/2];
sites=[];
for i=1:nx
    for j=1:ny
        v{i,j}=a1.*(i-1)+a2.*(j-1);
            sites=[sites,v{i,j}+a1/4,v{i,j}+a2/4,v{i,j}+(a2+a1)/4,v{i,j}-a1/4,v{i,j}-a2/4,v{i,j}-(a2+a1)/4];
    end
end

 figure(1)
 plot(sites(1,:),sites(2,:),'xr')
 hold on
 %%
 N_vertex=length(sites(1,:))
 Adj=zeros(N_vertex,N_vertex);
 diffs=zeros(N_vertex,N_vertex);
  Adj2=zeros(N_vertex,N_vertex);
  Adj3=zeros(N_vertex,N_vertex);



 for b1=-1:1
     for b2=-1:1
         sitesb{b1+2,b2+2}=sites+b1*nx*a1+b2*ny*a2;
     end
 end
 
for i=1:N_vertex
    for j=1:N_vertex
        for b1=1:3
            for b2=1:3
                ds2h(b1,b2)= abs(sites(1,i)-sitesb{b1,b2}(1,j))^2+abs(sites(2,i)-sitesb{b1,b2}(2,j))^2;
            end
        end
        ds2=min(ds2h(:));
        if sqrt(ds2)<0.5+0.01
            Adj(i,j)=1;
        elseif sqrt(ds2)<0.7
            Adj2(i,j)=1;
        elseif sqrt(ds2)<0.8
            Adj3(i,j)=1;
        end
        
         diffs(i,j)=sqrt(ds2);
    end
end


Adj=heaviside(Adj+Adj'-0.1);
Adj=Adj-diag(ones(1,length(Adj)));
for i=1:N_vertex
    for j=1:N_vertex
         diffs(i,j)=min(diffs(i,j),diffs(j,i));
    end
end
%Adj=Adj+Adj';
% hold on
%%
a=[0,1];
tic
for k=2:N_vertex
    tic
    k
    free=a;
    for j=1:k-1
        if Adj(k,j)==1
            free=free(bitget(free,j)==0);
        end
    end
    a=[a,free+2^(k-1)];
    tt(k)=toc
    dim(k)=size(a,2);
end
dimF=size(a,2)
toc

n=length(Adj)
toc
tic
s=zeros(size(a));
for k=1:n
s=s+bitget(a,k);
end
toc
%%
SX=cell(1,n);
PE=cell(1,n);
for k=1:n
[X,i1,i2]=intersect(a,bitxor(a,2^(k-1)));
SX{k}=sparse(i1,i2,ones(length(i1),1),dimF,dimF);
PE{k}=sparse(1:dimF,1:dimF,bitget(a,k),dimF,dimF);
end
HX=SX{1};
HD=PE{1};
for k=2:n
    HX=HX+SX{k};
    HD=HD+PE{k};
end
H2=0*PE{1};
H3=0*PE{1};
for k=1:n-1
    for j=k+1:n
        if Adj2(k,j)==1
        H2=H2+PE{k}*PE{j};
        end
        if Adj3(k,j)==1
        H3=H3+PE{k}*PE{j};
        end
    end
end

HR=0*PE{1};
for k=1:n-1
    for j=k+1:n
        if diffs(k,j)>0.5+0.01
        HR=HR+1/diffs(k,j)^6*PE{k}*PE{j};
        end
    end
end

%%
addpath('expmv')
%%
%tspan=0:0.01:10;
N1=300;
dt0=0.2;
s1=5*N1/100;
Delta0=-10;
Delta1=6;
Deltaspan=[linspace(0,Delta1-Delta0,N1),linspace(Delta1-Delta0,0,N1)];
Omegaspan=[ones(1,2*N1-5*ceil(s1)),zeros(1,5*ceil(s1))];
windowSize = ceil(20*N1/100);
b = (1/windowSize)*ones(1,windowSize);
Deltaspan=filter(b,1,Deltaspan)+Delta0;
Omegaspan=filter(b,1,filter(b,1,filter(b,1,Omegaspan)));
dtspan=dt0*ones(1,2*N1);
tvals=cumsum(dtspan);
psi=spalloc(dimF,1,1);
psi(1)=1;
for k=1:N1
    k
    dt=dtspan(k);
    H=-Deltaspan(k)*HD+Omegaspan(k)/2*HX;%+0*H2+0*0.66^6/0.75^6*H3+0.5*0*HR;%+0.2*H2;%+0.2*2^6/sqrt(3)^6*Hs3;
    psi=expv(dt,-1i*H,psi,10^(-6));
    [ev,ew]=eigs(H,50,'sa');
    EGS(:,k)=diag(ew);
    EV{k}=ev;
    densityGS(k)=ev(:,1)'*(HD)*ev(:,1);
    density(k)=psi'*(HD)*psi;
    pops(:,k)=abs(psi'*ev).^2;
end
%%
  figure(3)
  yyaxis left
  plot(cumsum(dtspan(1:N1)),Deltaspan(1:N1),'k')
  set(gca,'ycolor','k')
  hold on
  ylabel('$\Delta/\Omega_0$')
  yyaxis right
  plot(cumsum(dtspan(1:N1)),Omegaspan(1:N1),'-r')
  ylabel('$\Omega/\Omega_0$')
  ylim([0,1])
  xlim([0,N1*dt])
  hold off
  set(0,'defaulttextInterpreter','latex')
  xlabel('$\Omega_0 t$')
  set(gca,'FontSize',16)

%plot(Deltaspan, EGS-EGS(1,:))
%plot(Deltaspan./Omegaspan, (EGS(l,:)-EGS(1,:))./Omegaspan,'color', [pops(,0,0])
%hold on

%%
figure(5)
x = (Deltaspan(1:N1)./Omegaspan(1:N1));
y = ((EGS(:,1:N1)-EGS(1,1:N1))./Omegaspan(1:N1));
figure(5)
%plot(x,y,'-','linew',0.1,'color',[0.5,0.5,0.5])
hold on
  xlim([-6,6])
  ylim([0,5])

figure(5)
for l=length(ew):-1:1
    y = ((EGS(l,1:N1)-EGS(1,1:N1))./Omegaspan(1:N1));
    z = zeros(size(x));
col = log(pops(l,1:N1));  % This is the color, vary with x in this case.
surface([x;x],[y;y],[z;z],[col;col],...
        'facecol','no',...
        'edgecol','interp',...
        'linew',0.5);
  
    hold on

end
   
 colormap(flipud(cool))
 colormap(cool)
  xlim([-2,6])
  %ylim([-0.3,2.5])
  hold off
  caxis([-2,0])
  set(0,'defaulttextInterpreter','latex')
  xlabel('$\Delta/\Omega$')
  ylabel('$E_n-E_0$')
  set(gca,'FontSize',12)
  colorbar
  hold off
  %%
  figure(6)
  plot(x,densityGS./n,'k--')
  xlabel('$\Delta/\Omega$')
  ylabel('$\sum_i\langle n_i\rangle /N$')
  hold on
  plot(x,density./n,'r-')
  xlim([-4,4])
  xlabel('$\Delta/\Omega$')
    set(gca,'FontSize',18)
hold off
  figure(7)
  plot((x(1:end-1)+x(2:end))/2,diff(densityGS)./n,'k--')
    ylabel('$\frac{\partial}{\partial (\Delta/\Omega)}\sum_i\langle n_i\rangle /N$')
 xlabel('$\Delta/\Omega$')
 xticks([-5 -4 -3 -2 -1 0 1 2 3 4])
  hold on
  plot((x(1:end-1)+x(2:end))/2,diff(density)./n,'r-')
  xlim([-4,4])
    set(gca,'FontSize',18)
hold off
% ind0=find(hd==6);
% p0=zeros(1,dimF);
% p0(ind0)=1;
% ind1=find(hd==5);
% p1=zeros(1,dimF);
% p1(ind1)=1;
% sum(abs(p0'.*psi).^2)
% sum(abs(p1'.*psi).^2)
% sum(abs(p0'.*EV{end}(:,1)).^2)
% sum(abs(p1'.*EV{end}(:,1)).^2)
%%
indMIS=find(diag(HD)==max(HD(:)));
figure(9)
[val,ind]=sort(abs(EV{end}(indMIS,1)).^2,'descend');
plot(abs(EV{end}(indMIS(ind),1)).^2,'--xk')
xlabel('$i$ (MIS index)')
ylabel('$|c_i|^2$')
set(gca,'FontSize',14)
hold on
plot(abs(psi(indMIS(ind))).^2,'-xr')
figure(10)
plot(angle(psi(indMIS(ind))),'-xr')
ylim([-pi,pi])
ylabel('$\arg(c_i)$')
set(gca,'FontSize',14)